DETERMINANT APPROXIMATIONS 
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Abstract. A sequence of approximations for the determinant and its logarithm of a complex 
matrix is derived, along with relative error bounds. The determinant approximations are derived 
from expansions of det(X) = exp(trace(log(X))), and they apply to non-Hermitian matrices. Ex- 
amples illustrate that these determinant approximations are efficient for lattice simulations of finite 
temperature nuclear matter, and that they use significantly less space than Gaussian elimination. 

The first approximation in the sequence is a block diagonal approximation; it represents an 
extension of Fischer's and Hadamard's inequalities to non-Hermitian matrices. In the special case of 
Hermitian positive-definite matrices, block diagonal approximations can be competitive with sparse 
inverse approximations. At last, a different representation of sparse inverse approximations is given 
and it is shown that their accuracy increases as more matrix elements are included. 
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1. Introduction. For a complex matrix we present approximations for the de- 
terminant and its logarithm, together with error bounds. 

The approximations were motivated by a problem in computational quantum 
field theory: the simulation of finite temperature nuclear matter on a lattice [16] . In 
this application, the logarithm of the determinant is desired to 2-3 significant digits. 
The matrices are sparse, and non-Hermitian. Because the desired accuracy is low, 
LU decomposition with partial pivoting [12l §14.6], [2lJ §3.18] is too costly. Since 
the matrices are not Hermitian positive-definite, sparse approximate inverses }19j , 
Gaussian quadrature based methods [3] , and Monte Carlo methods [TH1 §4] or hybrid 
Monte Carlo methods [7J HO] do not apply. Monte Carlo and quadrature-based 
methods can be extended to non-Hermitian matrices, however then the sign of the 
determinant is usually lost, e.g. [2j §3.2.3]. 

To approximate the determinant det(M) we decompose M = Md + M Q g such 
that Mjj is a non-sing ular matrix. Then det(M) = det(M D ) det(I+M~~M off ), where 
/ is the identity matrix. In 

det(J + M D x M oS ) = exp(trace(log(7 + M^Mob))), 

we expand log(J + M g ), obtaining a sequence of increasingly accurate approx- 
imations. Error bounds for these approximations depend on the spectral radius of 
M^MoS. 

Overview. The determinant approximations and their error bounds are pre- 
sented in §S] Approximations from block diagonals f H2.ip are extended to a sequence 
of higher order approximations ( ^2.2[) . They simplify for checkerboard matrices f£ )2.3[) 
which occur in the neutron matter simulations in |16j . Comparisons with sparse in- 
verse approximations of determinants, which are limited to Hermitian positive-definite 
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matrices (jj3J) illustrate the competitiveness of block diagonal approximations. As 
expected, the accuracy of sparse inverse approximations increases as more matrix 
elements are included. Numerical results with matrices from nuclear matter simula- 
tions (jjU show that determinant approximations of desired accuracy can be obtained 
fast, in 1-3 iterations; and that they require significantly less space than Gaussian 
elimination (with partial or complete pivoting). 

Notation. The eigenvalues of a complex square matrix A are Xj(A), and the 
spectral radius is p(A) = maxj |Aj(A)|. The identity matrix is I, and A* is the conju- 
gate transpose of A. We denote by log(X) and exp(A) the logarithm and exponential 
function of a matrix X, and by ln(a;) and e x the natural logarithm and exponential 
function of a scalar x. 

2. Determinant Approximations. We present approximations to the deter- 
minant and its logarithm, as well as error bounds. 

2.1. Diagonal Approximations. We present relative error bounds for the ap- 
proximation of the determinant by the determinant of a block diagonal. Let M be a 
complex square matrix of order n partitioned k block matrix 
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where the diagonal blocks Mjj are square but not necessarily of the same dimension. 
Analogously, decompose M — Md + M Q ff into diagonal blocks Mo and off-diagonal 
blocks M ff, 
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(2.1) 



The block diagonal matrix Md is called a pinching of M [H §11.5]. We consider the 
approximation of det(M) by the determinant of a pinching, det (Md); and in particular 
bounds of the form det(M) < det(Mo). The matrices for which such bounds are 
known to hold are characterized by eigenvalue monotonicity of the following kind. 
A complex square matrix M is a r-matrix if [8l pp 156-57]: 

1. Each principal submatrix of M has at least one real eigenvalue. 

2. If Si is a principal submatrix of M and S\\ a principal submatrix of S\ then 
Amin(S'i) < A m i n (Sii), where A m i n denotes the smallest real eigenvalue. 

3. A min (M) > 0. 

The class of r-matrices includes Hermitian positive-definite, M-matrices and totally 
non- negative matrices [H pp 156-57], [TU Theorem 1]. 

Hadamard- Fischer Inequality. If M is a r-matrix then [8l Theorem 4.3] 



det(M) < det(Mo) 



(2.2) 



Strictly speaking, (|2.2p is called a Hadamard-Fischer inequality only for k = 2 (HI 
(0.5)]. If k = 2 and M is Hermitian positive-definite then (|2.2p is Fischer's inequal- 
ity [T31 Theorem 7.8.3]. If k = n and M Hermitian positive-definite then (|2.2p is 
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Hadamard's inequality [B] Theorem 8], JT3J Theorem 7.8.1]. Extensions of (|2.2[) to 
generalized Fan inequalities are derived in [T71 [T5] . The Hadamard-Fischer inequality 
(|2.2I) implies the obvious relative error bound for the determinant of a pinching, 

c det(Afo) - dct(A/) < 
det(Mr>) 

In the theorem below we tighten the upper bound. 

Theorem 2.1. Lef M be a complex matrix of order n. If det(M ) is real, Md 
is non-singular with det(Afo) real, and all eigenvalues Xj(M^ M Q g) are real with 
X j (M^ 1 M oS ) > -1, then 

det(M D ) - det(M) 



< , ' r , — - < 1 - e rrx— 

det(M D ) 

where p = p(M^ 1 M oS ) and A min = mini<j< n Xj(M^ M oS ). 

Proof Write det(M) = det(M D ) det(I + A), where A = M=; 1 M off . Since I + A 
is non-singular, [14j Theorem 6.4.15(a)] and [TH Problem 6.2.4] imply det(I + A) = 

exp(trace(log(7 + A))). Furthermore, log(7 + A) = J2™ =1 ' AP PS (7) in §9.8]. 

Hence 

det(J + A) = cxp M 1 + X A A )) 

Because Xj(A) > —1, 1 < j < n, we can apply the inequality < ln(l + A) < A [JJ 
4.1.33] to obtain 

exp(trace(A)) e 1 + A mi„ < dct(7 + A) < exp(trace(A)). 

At last use the fact that Md is block diagonal and trace(M Q gf) = trace(A) = 0. □ 

The upper bound for the relative error is small if the eigenvalues of M^ 1 M ff are 
small in magnitude but not too close to —1. The pinching det(Mo) can be a bad 
approximation to det(M) when / + M^Mog is close to singular. If det (Md) > 
then Theorem 1 2 . 1 1 implies a lower bound for det(M), 



.» det (M D ) < det(M) < det(M D )- 



In the argument of the exponential function in Theorem l2.1l we have A m i n < because 
M^Moff has a zero diagonal so that trace(M^ M Q s) = 0. Hence np 2 /(l + A min ) > 

9 

np . 

Corollary 2.2. Theorem \2.1\ holds for Hermitian positive- definite matrices. In 
particular, Theorem \2.1\ imvlies error bounds for Fischer's and Hadamard's inequali- 
ties. 

The following example shows that |det(M)| < |dct(M£>)| may not hold when 
Mp M a g has complex eigenvalues, or real eigenvalues smaller than — 1. 

Example 1. Even if all eigenvalues of M^ M Q g satisfy \Xj(Md M g)| < 1, it 
is still possible that | det(M)| > | det(M]o)| if some eigenvalues are complex. 

Consider 



1 a\ \ 0\ (0 a 
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Then Xj(M D 1 M a g) = ±a, and det(Af) = 1 — a 2 . Choose a — ^i, where i = 
\J —\. Then both eigenvalues of Mq M & are complex, and \\j(M D ~ 1 M g)\ < 1. But 
det(M) = 1.25 > 1 = det(Af£>). 

The situation det(Aft>) > det(M) can also occur when Mf^M s has a real eigen- 
value that is less than — 1 . If a — 3 in the matrices above then one eigenvalue of 
M^Matt is -2, and | dct(M)| = 8 > det(M D ) = 1. In general, | det(M)|/det(Mo) -> 
oo as \a\ — > oo. 

This example illustrates that, unless the eigenvalues of Mp M a g are real and 
greater than —1, det(Afo) is, in general, not a bound for det(M). In the case of 
complex eigenvalues, however, we can still determine how well det(Afo) approximates 
det(M). Below is a relative error bound for the case when M is 'diagonally dominant', 
in the sense that the eigenvalues of M^Mqs are small in magnitude. 

Theorem 2.3 (Complex Eigenvalues). Let M be a complex matrix of order n. 
If Md is non-singular and p = p(Mp M s) < 1 then 

|det(M)-det(M23)| 

-<cpe p , where c = — nm(l — p). 



|det(M D )| 
If also cp < 1 then 



|dct(M) - dct(Afp)| 7 
|det(Afo)| " 4 CP - 



Proof. This is a special case of Theorem 12.61 □ 

Corollary 2.4. Theorem \ 2.3\ holds for the following classes of matrices: M- 
matrices; Hermitian positive-definite matrices if k = n; Hermitian positive definite 
block tridiagonal matrices with equally- sized blocks of dimension n/k. 

Proof. In all cases p(M^ 1 M oS ) < 1. □ 

In the special case of strictly diagonally dominant matrices, Theorem 12 .31 leads to 
a bound for the approximation of det(Af) by the product of the diagonal elements. 

Corollary 2.5. If the complex square matrix M = (ra, 3 )i<, j< n is strictly row 
diagonal dominant then 



' n I < cpe p , where p < max > 

lli-l m ii 1 



-nln(l — p). 



If also cp < 1 </ien 



in"=i 
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Proof. This is a consequence of Gerschgorin's theorem [H Theorem 7.2.1]. □ 
Corollary 12 . 5 1 implies that the product of diagonal elements is a good approxima- 
tion for det(M) if M is strongly diagonally dominant. 

2.2. A Sequence of General Higher Order Approximations. We extend 
the diagonal approximations in §2.1l to a sequence of more general approximations that 
become increasingly more accurate. These approximations, called 'zone determinant 
approximations' in |16j . are justified in the context of nuclear matter simulations. As 
before, decompose M = Md + M g into diagonal blocks Md and off-diagonal blocks 
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M g (actually, our results hold for any decomposition M — M + Me where M is 
non-singular and p{M^ 1 Me) < 1). Below we give a sequence of approximations S m 
for ln(det(M)) and A m for det(M), as well as absolute bounds for S m and relative 
bounds for A m . An absolute bound for the logarithm suffices because ln(det(M)) > 1 
in our applications. 

Theorem 2.6. Let M = Md + M Q s be a complex matrix of order n, Md be 
non-singular and p = p{M D \ l M Q g) < 1. Define 

S m = ln(det(AfD)) + ~ — " trace((Mo 1 M off )P), A m = e 5 '", m > 1. 

Then 

| ln(det(M)) - S m \ < cp m , I det(M) - A m | < ^ m ^ 

\^m | 

where c = — nln(l — p). If also cp 771 < 1 then 

|det(M)-A m | 7 m 
|A m | "4 CP ' 



Proo/. As in the proof of Theorem |2~T1 det(M) = det(M D ) det(I + A), where 

oo (-l)P 

p=l p 



A = M^Moff and log(7 + A) = E^li ^^ p - Hence 



00 f— 

trace (log(7 + A)) = V ^ — '- trace(4 p ). 

1 P 

P =i 



Define the truncated sums 



L m ee V ^ i trace(A p ), D m = e L "\ m > 1. 



Then 



trace(log(7 + A)) =L m + z, z = £ \ ln(l + \{A)) - ]T \{Ay 



i=l k P=l ^ 



Applying to each of the n terms the inequality 

<-ln(l-|A|)|Ar 



ln(l + A) - £ A* 

P =i 



[2 4.1.24], [TJ 4.1.38] gives | ln(det(J + A)) - L m \ < cp m . The first bound follows now 
with S m = ln(det(Afo)) + L m . 

From the first bound, the fact that det(J + A) = D m e z , and \e z - 1| < |z| e |z| [TJ 
4.2.39] follows 

\det(i + A)-ru <cpm eCp ™ 



\D m \ 

We get the second bound from A m = det(Mo) D r , 
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If also cp m < 1 then Q] 4.2.38; 



6eb{I + A)-D m \ < 7_ cpn 



I An | " 4 

□ 

The accuracy of the approximations in Theorem l2.6l is determined by the spectral 
radius p of M^Mog. In particular, the absolute error bound for the approximation 
d m is proportional to p m , hence the approximations tend to improve with increasing 
m. The numerical results in Sections [3] and |4] illustrate that the pessimistic factor in 
the bound | ln(det(M)) — S m \ < — nln(l — p) p m is n. We found that replacing n by 
the number of eigenvalues whose magnitude is close to p makes the bound tight. The 
approximations for the logarithm can be determined from successive updates 

f— 1 "> 1,1-1 

S = ln(det(Mr>)), S m = S m -i + - — trace^M" 1 M oS ) m ), m > 1, 

m 

and A m = e 5m . Note that e 5 ° = det(Mo) is the block diagonal approximation from 
(|2.ip . Hence Theorem 12.31 is a special case of Theorem 12.61 with m = 1 . If a block 
diagonal determinant approximation is sufficiently accurate, as in <21 it can be much 
cheaper to compute than a determinant via Gaussian elimination. 

2.3. Checkerboard Matrices. For this particular class of matrices, which oc- 
curs in our applications [TB] , every other determinant approximation A m has increased 
accuracy. We call a matrix M with equally sized blocks Mij of dimension n/k in (12 . 1 [) 
an odd checkerboard matrix (with regard to the block size n/k) if Mij = for i and j 
both even or both odd, 1 < i, j < k; and an even checkerboard matrix if Mij = for i 
odd and j even or vice versa. An odd checkerboard matrix has zero diagonal blocks, 
hence its trace is zero. 

Theorem 2.7. //, in addition to the conditions of Theorem \2.b\ M g is an odd 
checkerboard matrix then 



So = ln(det(M D )), 




if m is odd 

, ( (M- 1 M oii ) m \ . t 

trace - — — if m is even. 

V m J 



Proof. If A and B are odd checkerboard matrices (with regard to the same block 
size) then AB is an even checkerboard matrix. If A is an odd checkerboard matrix 
and B an even checkerboard matrix then AB and BA are odd checkerboard matrices. 

Since M^ 1 M Q ff is an odd checkerboard matrix, so are the powers (M^ 1 M a g) p for 
odd p. This means trace ((M^ 1 M a g) p ) = for odd p. Hence the approximations in 
Theorem 12.61 satisfy 5 m = 5 m -i for m odd. For m even 

□ 

Theorem 12.71 shows that an odd-order approximation is equal to the previous 
even-order approximation. Hence the even-order approximations gain one order of 
accuracy. 
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3. Comparison with Sparse Inverse Approximations. In the special case of 
Hermitian positive-definite matrices, we illustrate that block-diagonal determinant ap- 
proximations (see Corollary 12. 2p can compare favourably with approximations based 
on sparse approximate inverses |I9j . We also show that the accuracy of sparse inverse 
approximations increases when more matrix elements are included. 

Idea. To understand how sparse inverse approximations work, we first consider a 
representation of the determinant based on minors of the inverse p~3l §0.8.4]. If M is 
Hermitian positive-definite of order n, and is the leading principal submatrix of 
order i of M, then 13 §0.8.4] det(M) = det(M„_i)/cr„, where a n = (M _1 ) n „ is the 
trailing diagonal element of M _1 . Using this expression recursively for det(M n _i) 
gives 

n 1 

det(M) = J| — , where en = (A^ 1 ),,. 
i=i Ut 

Determinant approximations based on sparse approximate inverses replace leading 
principal submatrices Mi by just principal submatrices Si. Specifically [TU §3.2], let 
M be Hermitian positive-definite, and let Si be a principal submatrix of Mi, such that 
Si includes at least row i and column i of M. The two extreme cases are Si — ma 
and Si = Mi. In any case, ma is the trailing diagonal element of Si, i.e. S nuni — ma, 
where rn is the order of Si, 1 < rii < i. Let o~i be the trailing diagonal element of 
S^~ , i.e. (Tj = (S^ 1 )n i ,n i - In particular o\ = . Given n such submatrices Si, 
1 < i < n, the sparse inverse approximation of det(M) is defined as [TH1 Algorithm 
3.3]. 

n -. 

a = H-. (3.1) 

i— i 

The sparse approximate inverse method performs Cholesky decompositions Si — 
LiL*, where L; is lower triangular, 2 < i < n, and computes l/<7j = ((Li) nuni ) 2 . 

Monotonicity. We show monotonicity of the sparse inverse approximations in 
the following sense: If the dimensions of the submatrices Si are increased then the 
determinant approximations can only become better. 

Lemma 3.1. If 

m k 

M =k U* s) 

is Hermitian positive- definite then (S^ 1 )u < (Af _ ) m +i,m+i) 1 < « < fc- 

Proof. The proof follows from [51 (4)] and the Shermann-Morrison formula [HI 
(2.1.4)]. □ 

Lemma 13.11 implies the following lower and upper bounds for sparse inverse ap- 
proximations; the lower bound was already derived in |19l (3.25)]. 

Corollary 3.2. If M is Hermitian positive-definite and a is a sparse inverse 
approximation in VS. 1\) then 

det(M) < a < Y[m li . 
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Corollary 13.21 implies that the product of diagonal elements cannot approximate 
the determinant more accurately than a sparse inverse approximation. Another con- 
sequence of Lemma |3. II is the monotonicity of the sparse inverse approximation in the 
following sense: If a principal submatrix Sj is replaced by a larger principal submatrix 
Sj then the determinant approximation can only become better. 

Theorem 3.3. Let M be Hermitian positive-definite of order n. If for some 
1 < j < n, Sj is a principal submatrix of Mj, and in turn Sj is a principal submatrix 
of Sj then 



det ( M)<ni<± n i 



where <jj is the trailing diagonal element of Sj . 

The next example of block diagonal matrices illustrates that sparse inverse ap- 
proximations can be inaccurate, even when sparsity is exploited to full extent. 

Block- Diagonal Matrices. Let 

(T 3 \ 

M = 

V t 3 ) 

be a block diagonal matrix of order n — 3k with n/3 diagonal blocks 

'3/2 -1 
T 3 = | -1 3/2 



The obvious block diagonal approximation (|2.1 



minant det(Mo) = det(M) = det(T 3 )"/ 3 
mation (|3.1[) we choose the submatrices 




'(*-l)(n/3)+l 



3/2, 



S t 



(»-l)(n/3)+2 



= n/3 gives the exact deter- 
(3/8)™/ 3 . For the sparse inverse approxi- 

3/2 



-1 



3/2 



S 



i(n/3)j 



1 < i < n/3. 



The sparse inverse approximation of det(Ts) is det(T3) + 2/3. It has no accurate digit 
because the relative error is 16/9. The sparse inverse approximation of det(M) is 
a = (det(T 3 ) + 2/3))™ /3 . For instance, when n = 300 then det(M) 4 • 10 17 while 
the sparse inverse approximation gives a sa 4 • 10 33 . 

Tridiagonal Toeplitz Matrices. A block diagonal approximation can be more ac- 
curate than a sparse inverse approximation if the dimension of the blocks is larger 
than 1. Let 

/ 2 -1 \ 



Tn 



-1 



V 



-1 



-1 

2/ 



be of order n; then det(T„) = n + 1. In the sparse inverse approximation (|3.1|i 
we fully exploit sparsity by choosing S\ = 2 and Si = T%, 2 < i < n; hence the 
approximation is a = 2 (3/2)" -1 . When Mb in (|2.1[) consists of k equally sized blocks 

of dimension n/k then det(Mu) = (det(T„/fc)) fc = ((n/k) + l) k . For a block size 
n/k > 4, det(Mo) < it, so the block diagonal approximation is more accurate than 
the sparse inverse approximation. 



2-D Laplacian. We show that for this matrix both, the block-diagonal and the 
sparse inverse approximations are accurate to at most one digit. The coefficient matrix 
from the centered finite difference discretization of Poisson's equation is a Hermitian 
positive-definite block tridiagonal matrix 111] 9.1.1] 



M 



V 



-In 

T 

± in 



\ 



-Ir. 

T m ) 



where T„ 



/4 -1 
-1 4 

V 



-1 

4/ 



Here T m is of order m, and M is of order n = m 2 (note that the matrix considered in 
[19} §5] equals (n 4- 1) 2 M). The exact determinant is [11, Theorem 9.1.2] 



det(M) 




sin 



2(m+ 1) 



sin 



2(m+ 1), 



We compute the logarithm of this expression and compare it to the approximations. 
A block diagonal approximation (|2.1[) with k = m gives det(Mo) = det(T m ) m , where 
[H §28.5] 



1 



det(T m ) = TT ( 4 + 2cos- 

»=1 v 

If the matrices in the sparse inverse approximation p. II) are 

-r 



5i =4, 



5,; = 



4 
-1 



2 < i < m+1, 



Sj = 





4 

-1 



m+2 < j < n, 



then l/<7i = 4, l/o* = 15/4, 2 < i < m + 1, and l/cr, = 7/2, m + 2 < j < n. Thus 
the sparse inverse approximation is a = 4(15/4) m (7/2) Tl ~ m ~ 1 . 

Table EH] lists errors for the block diagonal and sparse inverse approximations for 
n = 900, n = 10000 and n = 40000. Columns 3 and 4 represent the relative errors 

| ln(det(M ))-ln(det(M))|/| ln(det(M))| and | ln(tr)-La(det(M))|/| ln(det(M))|, 

while columns 5 and 6 represent the relative errors 

|dct(M D ) 1 /"-det(Af) 1 /"|/|det(Af) 1 /™| and |<r 1 / T, -det(Af) 1 / n |/| det(M) 1 / n |. 

We include the last two errors to allow a comparison with the approximation of 
det ((n + l) 2 M) 1/n in Table 5.1]. The table shows that all relative errors lie 
between 0.06 and 0.2. Hence both approximations, block diagonal and sparse inverse, 
are accurate to at most one significant digit. To estimate the tightness of the bound 

| ln(dct(M D )) - In det (M) | < (-n | ln(l - p)\) p 

in Theorem l2.6[ consider the case n = 900. Here p{M^-M oS ) « .9898 and | ln(l— p)| « 
4.5845. The true error is 

| ln(det(MD)) - lndet(M)| w 122.4966 w 261n(l - p) p. 

The matrix M~^M a s has 26 eigenvalues with magnitude at least 0.9. Thus the 
pessimism of the bound comes from the factor n. 
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n 


ln(det(M)) 


rel. error 
in hx(M D ) 


rel. error 
in ln(cr) 


rel error 

ji.fl/n 

m M D 


rel. error 
in cr 1 /" 


900 
10000 
40000 


1.0650e+03 
1.1717e+04 
4.6761c+04 


0.1150 
0.1246 
0.1269 


0.0607 
0.0698 
0.0719 


0.1458 
0.1572 
0.1599 


0.0745 
0.0852 
0.0877 



Table 3.1 

Errors in the block diagonal approximation Mo and the sparse inverse approximation a for the 
Laplacian. 



4. Application to Neutron Matter Simulations. In [IB] we consider the 
quantum simulation of nuclear matter on a lattice, and in particular how to calcu- 
late the contribution of nucleon-nucleon-hole loops at non-zero nucleon density. The 
resulting method, called zone determinant expansion, is based on the sequence of 
approximations in Theorem 12.61 Here we illustrate that 3 iterations of the zone de- 
terminant expansion give an approximation accurate to 3 digits, and that the method 
uses less space than a determinant computation based on Gaussian elimination (with 
partial or complete pivoting). 

In [16] we derive a particle interaction matrix M whose determinant det(M) is 
not positive, and complex in general. Hence stochastic methods such as hybrid Monte 
Carlo methods [3 [TOj [20] do not give the correct sign or phase of det(M). This was 
the motivation for approximating ln(dct(M)) via a zone determinant expansion, i.e. 
Theorem 12.61 Below we discuss the structure of M and a physically appropriate zone 
determinant expansion. 

The particle interactions are considered on a 4-dimensional lattice (3 dimensions 
for space and one for time). Let the dimensions of the lattice be L x L x L x L t , where 
L t represents the time direction. Also let the number of particles per lattice point be 
s. Then the interaction matrix M has dimension nxn where n = L 3 L t s. We partition 
the lattice into separate spatial zones (or cubes) of dimension m x m x m (constraints 
on m are discussed in [16]). Therefore particle interactions between any two zones are 
represented by matrix blocks of dimension m 3 L t s. As a consequence, it makes sense 
to approximate det(M) by the product of principal minors associated with particle 
interactions inside spatial zones. Without loss of generality we assume that the lattice 
points are ordered such that the submatrix Mij of order m 3 L t s represents particle 
interactions between zones i and j. With k = (L/m) 3 this gives the partitioning 
M — Md + M s in (|2.1j) . where Mo represents particle interactions in the zone 
interiors, while M g represents interactions among different zones. In [16) we explain 
that the spectral radius p = p(M D ~ 1 M ff) can be reduced by increasing the dimension 
m of the spatial zones. 

We illustrate the zone expansion on a small lattice simulation, where we can 
compare the approximations to the exact determinant. Specifically we consider the 
interactions between neutrons and neutral pions, on a 4 3 x 4 grid. The order of the 
interaction matrix M is 4 3 x 4 x 2 = 512. Its properties are listed in Table |4~T1 

In the context of the particular application in [16] . we can partition the lattice 
into zones with dimension m = 1. The resulting partitioning has blocks My = 
M 8 ( i _ 1 < )+1 . 8L8 (j_ 1 * )+ i.j, 1 < i,j < 64, of dimension 4x2 = 8. Thus k — 64 in the 
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nz = 4608 



Fig. 4.1. Sparsity structure of the interaction matrix M . 



order 


n = 512 




number of non-zeros 


9n 


see Figure 14.11 


structure 


complex non-Hermitian 


see Figure 14.21 


norm 


||M|| F « 49.5 




condition number 


IIMHxIlA/- 1 !!! « 177 


condcst(-) in MATLAB 6 


non-normality 


\\M*M - MM*\\ F rs 57 




eigenvalues 


complex 


see Figure [4731 


determinant 


det(M) = 8.5361 • 10 65 + 1.4 1 68 • 10 64 i 
ln(det(M)) = 151.81 + 0.016599« 


det(-) in Mat lab 6 



Table 4.1 
Properties of the interaction matrix M . 



10 20 30 40 50 60 J 

nz = 448 ' ' ' ' n ,.;, 5 " ' a s 

Fig. 4.2. Non-zero 8x8 blocks in the interaction matrix M, and sparsity structure of a single 
8x8 diagonal block. 
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x + a. +i +3* 

A + /+ if 
* ++* + ++ + 
+ + + ^ + 

+*++ + *+ 
+ 1 + ++ 



#4 +*- 



+ ++*++ 



-^lHl^iljl|iiii^ 



Fig. 4.3. Eigenvalue distribution of the interaction matrix M. 





Fig. 4.4. Sparsity structure of the matrices M D ^M f[ and (M^M^) 2 



block diagonal approximation (|2.ip . Figure l4~2l shows the distribution of the 448 blocks 
with non-zero elements. Each diagonal block Ma contains 24 non-zero elements, its 
sparsity structure is shown in Figure 14.21 

The zone partitioning is bipartite, i.e. My = for i and j both even or both 
odd, and i ^ j, 1 < i,j ' < fc. Therefore M fi is an odd checkerboard matrix. Figure 
4.51 illustrates this checkerboard pattern in the leading principal submatrix of order 
32 of Mp Mag. The sparsity structures of the matrices Mp M a g and (M^ 1 M g) 2 
is shown in Figure 14.41 Because of the checkerboard structure Theorem 12.71 implies 
trace((M i ^ 1 M ff) p ) = for odd p, and 5 p -x = 5 p . Table l4~2l therefore contains only 
approximations of even order. 

Table 14.21 shows errors in the approximations 5j and Aj for approximations up to 
order 8. Columns 2, 3 and 4 represent the absolute errors 



|»(ln(det(M))) - $t(5j)\, |3f(ln(det(M))) - 3(^)1, 
Columns 6 and 7 represent the relative errors 



ln(det(M))-^|. 



|ln(det(M))-^|/|^| 

The spectral radius p = p(M^ 1 M Q g) 
of Theorems O and is c « 554. 



and |det(M)-Aj|/|Aj|. 
.6613, and the constant in the error bounds 
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Fig. 4.5. Sparsity structure of the leading principal submatrix of order 32 of M D ^ M a g. 



3 


abs. error 


abs. error 


abs. error 


P> 


rel. error 


rel. error 




in ft(6j) 


in ^s(Sj) 


in Sj 




in Sj 


in Aj 





5.1000 


0.0017 


5.1000 




0.0348 


163.0282 


2 


0.4817 


0.0025 


0.4817 


0.4374 


0.0032 


0.3823 


4 


0.0909 


0.0016 


0.0909 


0.1913 


0.0006 


0.0951 


6 


0.0225 


0.0008 


0.0226 


0.0837 


0.0001 


0.0223 


8 


0.00665 


0.0003 


0.0066 


0.0366 


0.00004 


0.0066 



Table 4.2 

Errors in the approximations Sj and Aj for the interaction matrix M . 



Table l4~2l illustrates that |ln(det(M)) — Sj\ ps p> , i.e. the absolute errors in the 
logarithm are almost proportional to the powers of the spectral radius of M~^M a s ■ In 
this case the constant c is too pessimistic, because many eigenvalues of M z ^ 1 M g have 
magnitude much less than p. For instance, 160 eigenvalues of M^ 1 M r have mag- 
nitude 10 -15 . The imaginary parts of the logarithms appear to converge faster than 
the real parts. The block diagonal approximation So = ln(det(M£>)) for ln(det(M)) 
has an accuracy of 2 digits. Two more iterations give an approximation S2 that is 
accurate to 3 digits. 

We briefly compare the computation of <5o and S2 to a determinant computation 
by Gaussian elimination of M. Gaussian elimination with complete pivoting gives 
PMQ = LU , where P and Q are permutation matrices, L is unit lower triangular 
and U is upper triangular. Figure I4.6[ which shows the sparsity structure of the 
matrices PMQ, L and U, illustrates that Gaussian elimination with complete pivoting 
completely destroys the sparsity structure of M. The matrices L and U together have 
about 162n non-zeros, compared to 9n in M. In contrast, the determinant expansion 
requires no significant additional space for 6q; and 48n non-zeros for M~^M a s and 
n non-zeros for the trace of (Mp M Q s) 2 . That's (48 + l)n = 49n non-zeros, about 
one third of the non-zeros produced by Gaussian elimination with complete pivoting. 
Gaussian elimination with partial pivoting essentially preserves the sparsity structure 
of M but produces 342n non-zeros. 
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